
!-------------------------------------------------------------------------------------
subroutine read_input()
!-------------------------------------------------------------------------------------

  use dimensions
  use mpi_wrap
  use parameters


  implicit none

  logical :: filefound
  integer :: fin, fout

  real    :: u, uAngle, theta
  integer :: iDir, i, imax

  fin  = 1
  fout = 2


 !-- print header to output file --

  write(msg,*) char(10), ' Reading input file'
  call mpi_wrap_msg



  !-- check if input file exist --

  inquire(file = trim(runname)//'.in', exist=filefound)
  if(.not.filefound) then
     write(msg,*)    'cannot find the input file'
     call mpi_wrap_err
  end if

  open(fin,file = trim(runname)//'.in',form='formatted')

  !-- ignore header --

  read(fin,*) 
  read(fin,*) 
  read(fin,*)
  read(fin,*) 
  read(fin,*)

  !-- background electric field --

  read(fin,*) 
  read(fin,*) 
  read(fin,*)

  read(fin,*,ERR=502,END=502)  phi0
  read(fin,*,ERR=502,END=502)  kz
  read(fin,*,ERR=502,END=502)  vPhase
  read(fin,*,ERR=502,END=502)  SpeckleWidth
  read(fin,*,ERR=502,END=502)  SpeckleDepth
  read(fin,*,ERR=502,END=502)  SideLoss
  read(fin,*,ERR=502,END=502)  CollisionRate
  read(fin,*,ERR=502,END=502)  iE0Modulation


  !-- domain size --

  read(fin,*) 
  read(fin,*) 
  read(fin,*)

  read(fin,*,ERR=502,END=502)  nz
  read(fin,*,ERR=502,END=502)  Lx
  read(fin,*,ERR=502,END=502)  Ly
  read(fin,*,ERR=502,END=502)  Vmax


  !-- initial conditions --

  read(fin,*) 
  read(fin,*) 
  read(fin,*)

  read(fin,*,ERR=502,END=502)  u
  read(fin,*,ERR=502,END=502)  uAngle
  read(fin,*,ERR=502,END=502)  iModx
  read(fin,*,ERR=502,END=502)  iMody
  read(fin,*,ERR=502,END=502)  iModz
  read(fin,*,ERR=502,END=502)  ModAmp


  !-- numerical parameters --

  read(fin,*) 
  read(fin,*) 
  read(fin,*)

  read(fin,*,ERR=502,END=502)  restart
  read(fin,*,ERR=502,END=502)  nt
  read(fin,*,ERR=502,END=502)  dt
  read(fin,*,ERR=502,END=502)  iDiag
  read(fin,*,ERR=502,END=502)  iData
  read(fin,*,ERR=502,END=502)  Dxxxx
  read(fin,*,ERR=502,END=502)  Dzzzz
  read(fin,*,ERR=502,END=502)  Dvvvv
  read(fin,*,ERR=502,END=502)  Vcutoff
  read(fin,*,ERR=502,END=502)  kFilterAngle
  read(fin,*,ERR=502,END=502)  iuModulation
  read(fin,*,ERR=502,END=502)  iTransAdvect
  read(fin,*,ERR=502,END=502)  iLinHydro
  read(fin,*,ERR=502,END=502)  iLinVlasov
  read(fin,*,ERR=502,END=502)  iBcPerpPeriodic
  read(fin,*,ERR=502,END=502)  iEraseNonlin
  read(fin,*,ERR=502,END=502)  modeSave  
  read(fin,*,ERR=502,END=502)  iWaveFrame
  read(fin,*,ERR=502,END=502)  iFourier
  read(fin,*,ERR=502,END=502)  iDealias

  close(fin)


!-------------------------------------------------------------------------------------
! set pi and backgound velocities
!-------------------------------------------------------------------------------------

  pi = acos(-1.0)

  imax = jDir*npDir

  do iDir = 1, jDir
     i = myid * jDir + iDir - 1
     theta = 2*pi*i/imax + pi*uAngle
     ux0(iDir) = u * cos(theta)
     uy0(iDir) = u * sin(theta)
  enddo


!-------------------------------------------------------------------------------------
!  check input
!-------------------------------------------------------------------------------------

  if (myid .ne. 0) return

  open(fout, file=trim(runname)//'.out',  position='append')

  write(fout, *)
  write(fout, *)"Precompiled system size:"
  write(fout, *)
  write(fout, "('    jx =              ', (I6) )" ) jx
  write(fout, "('    jy =              ', (I6) )" ) jy
  write(fout, "('    jz =              ', (I6) )" ) jz
  write(fout, "('    jv =              ', (I6) )" ) jv
  write(fout, "('    jDir =            ', (I6) )" ) jDir*npDir

  write(fout, *)
  write(fout, *)"Input parameters:"
  write(fout, *)
  write(fout, "('    phi0 =            ', (e19.8) )" ) phi0
  write(fout, "('    kz =              ', (e19.8) )" ) kz
  write(fout, "('    vPhave =          ', (e19.8) )" ) vPhase
  write(fout, "('    SpeckleWidth =    ', (e19.8) )" ) SpeckleWidth
  write(fout, "('    SpeckleDepth =    ', (e19.8) )" ) SpeckleDepth
  write(fout, "('    SideLoss =        ', (e19.8) )" ) SideLoss
  write(fout, "('    CollisionRate =   ', (e19.8) )" ) CollisionRate
  write(fout, "('    iE0Modulation =       ', (L) )" ) iE0Modulation

  write(fout, *)
  write(fout, "('    nz =              ', (I6) )" )    nz
  write(fout, "('    Lx =              ', (e19.8) )" ) Lx
  write(fout, "('    Ly =              ', (e19.8) )" ) Ly
  write(fout, "('    Vmax =            ', (e19.8) )" ) Vmax

  write(fout, *)
  write(fout, "('    u =               ', (e19.8) )" ) u
  write(fout, "('    uAngle =          ', (e19.8) )" ) uAngle
  write(fout, "('    iModx =               ', (L) )" ) iModx
  write(fout, "('    iMody =               ', (L) )" ) iMody
  write(fout, "('    iModz =               ', (L) )" ) iModz
  write(fout, "('    ModAmp =          ', (e19.8) )" ) ModAmp

  write(fout, *)
  write(fout, "('    restart =             ', (L) )" ) restart
  write(fout, "('    nt =              ', (I6) )" )    nt
  write(fout, "('    dt =              ', (e19.8) )" ) dt
  write(fout, "('    iDiag =           ', (I6) )" )    iDiag
  write(fout, "('    iData =           ', (I6) )" )    iData
  write(fout, "('    Dxxxx =           ', (e19.8) )" ) Dxxxx
  write(fout, "('    Dzzzz =           ', (e19.8) )" ) Dzzzz
  write(fout, "('    Dvvvv =           ', (e19.8) )" ) Dvvvv
  write(fout, "('    Vcutoff =         ', (e19.8) )" ) Vcutoff
  write(fout, "('    kFilterAngle =    ', (e19.8) )" ) kFilterAngle
  write(fout, "('    iuModulation =        ', (L) )" ) iuModulation
  write(fout, "('    iTransAdvect =        ', (L) )" ) iTransAdvect
  write(fout, "('    iLinHydro =           ', (L) )" ) iLinHydro
  write(fout, "('    iLinVlasov =          ', (L) )" ) iLinVlasov
  write(fout, "('    iBcPerpPeriodic =     ', (L) )" ) iBcPerpPeriodic
  write(fout, "('    iEraseNonlin =        ', (L) )" ) iEraseNonlin
  write(fout, "('    modeSave =        ', (I6) )" )    modeSave
  write(fout, "('    iWaveFrame =      ', (I6) )" )    iWaveFrame
  write(fout, "('    iFourier =        ', (I6) )" )    iFourier
  write(fout, "('    iDealias =            ', (L) )" ) iDealias

  write(fout, *)

  close(fout)


  if(   (iFourier.ne.1).and.(iFourier.ne.-1).and.(iFourier.ne.0)  ) then
     write(msg,*) 'input variable iFourier must be in (-1, 0, 1)'
     call mpi_wrap_err
  endif

  if( nz > jz/2 -1) then
     write(msg,*) 'nz is too big'
     call mpi_wrap_err
  end if

  return

502 continue

  write(msg, *) 'cannot process input file'
  call mpi_wrap_err

end subroutine read_input

!-------------------------------------------------------------------------------------
